﻿using System;
using System.Collections.Generic;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace SSIM
{
    public class SSIMCalculator
    {
        //滑动窗口算法 
        private const int WindowSize = 11;//滑动窗口大小包含11个元素
        private const double GaussianSigma = 1.5;
        private readonly double[,] _gaussianKernel; //二维数组 高斯核数学结构

        private const double C1 = (0.01 * 255) * (0.01 * 255);
        private const double C2 = (0.03 * 225) * (0.03 * 225);

        public SSIMCalculator()
        {
            _gaussianKernel = GenerateGaussianKernel(WindowSize, GaussianSigma);
        }

        /// <summary>
        /// 生成一个size * size 的高斯核
        /// </summary>
        /// <param name="size">核尺寸必须是奇数，如11*11</param>
        /// <param name="sigma">高斯分布的标准差，控制模糊强度</param>
        /// <returns></returns>
        private double[,] GenerateGaussianKernel(int size, double sigma)
        {
            var kernel = new double[size, size];
            double sum = 0.0;
            int radius = size / 2; //核的半径

            //高斯公式，根据像素距离中心的偏移量(i,j)计算权重
            for (int i = -radius; i <= radius; i++)
            {
                for (int j = -radius; j <= radius; j++) 
                {
                    double value = Math.Exp(-(i * i + j * j) / (2 * sigma * sigma));
                    kernel[i + radius, j + radius] = value;
                    sum += value;
                }
            }

            // 归一化 
            //目的：确保高斯核的总数和为1，避免模糊操作改变图像亮度
            //操作：将所有权重除以sum，使权重分布满足概率性质
            for (int i = 0; i < size; i++)
            {
                for (int j = 0; j < size; j++)
                {
                    kernel[i, j] /= sum;
                }
            }

            return kernel;
        }


        // 将图像转换为灰度数组
        private byte[] ConvertToGrayscale(Bitmap image)
        {
            byte[] gray = new byte[image.Width * image.Height];
            int index = 0;

            for (int y = 0; y < image.Height; y++)
            {
                for (int x = 0; x < image.Width; x++)
                {
                    Color pixel = image.GetPixel(x, y);
                    byte value = (byte)(0.299 * pixel.R + 0.587 * pixel.G + 0.114 * pixel.B);
                    gray[index++] = value;
                }
            }

            return gray;
        }

        // 计算SSIM
        public double CalculateSSIM(Bitmap img1, Bitmap img2)
        {
            if (img1.Width != img2.Width || img1.Height != img2.Height)
                throw new ArgumentException("图像尺寸不匹配");

            byte[] gray1 = ConvertToGrayscale(img1);
            byte[] gray2 = ConvertToGrayscale(img2);

            int width = img1.Width;
            int height = img1.Height;
            int radius = WindowSize / 2;
            double totalSSIM = 0.0;
            int validPixels = 0;

            // 遍历每个像素的窗口
            for (int y = radius; y < height - radius; y++)
            {
                for (int x = radius; x < width - radius; x++)
                {
                    double mu1 = 0.0, mu2 = 0.0;
                    double sigma1Sq = 0.0, sigma2Sq = 0.0, sigma12 = 0.0;

                    // 计算窗口内的加权均值、方差和协方差
                    for (int i = -radius; i <= radius; i++)
                    {
                        for (int j = -radius; j <= radius; j++)
                        {
                            double weight = _gaussianKernel[i + radius, j + radius];
                            int index = (y + i) * width + (x + j);
                            double pixel1 = gray1[index];
                            double pixel2 = gray2[index];

                            mu1 += pixel1 * weight;
                            mu2 += pixel2 * weight;
                        }
                    }

                    for (int i = -radius; i <= radius; i++)
                    {
                        for (int j = -radius; j <= radius; j++)
                        {
                            double weight = _gaussianKernel[i + radius, j + radius];
                            int index = (y + i) * width + (x + j);
                            double pixel1 = gray1[index];
                            double pixel2 = gray2[index];

                            sigma1Sq += weight * (pixel1 - mu1) * (pixel1 - mu1);
                            sigma2Sq += weight * (pixel2 - mu2) * (pixel2 - mu2);
                            sigma12 += weight * (pixel1 - mu1) * (pixel2 - mu2);
                        }
                    }

                    // 计算当前窗口的SSIM
                    double numerator = (2 * mu1 * mu2 + C1) * (2 * sigma12 + C2);
                    double denominator = (mu1 * mu1 + mu2 * mu2 + C1) * (sigma1Sq + sigma2Sq + C2);
                    double ssim = numerator / denominator;

                    totalSSIM += ssim;
                    validPixels++;
                }
            }

            return totalSSIM / validPixels;
        }


    }
}
